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Abstract 



We propose an approximation to full relativity that captures the main grav- 
itational effects of dynamical importance in supernovae. The conceptual link 
between this formalism and the Newtonian limit is such that it could likely be 
implemented relatively easily in existing multidimensional Newtonian gravi- 
tational hydrodynamics codes employing a Poisson solver. As a test of the 
formalism's utility, we display results for rapidly rotating (and therefore highly 
deformed) neutron stars. 
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I. INTRODUCTION 



The collapsed cores of massive stars are relativistic bodies. Order-of-magnitude esti- 
mates of the gravitational potential, infall and outflow velocities, equatorial linear velocity, 
and microscopic nucleon velocity — computed using typical values of mass, radius, and ro- 
tation periods of pulsars — suggest that relativity cannot be neglected. More importantly, 
detailed comparisons of Newtonian and relativistic supernova simulations in spherical sym- 
metry show more compact cores and higher neutrino luminosities and average energies in 
relativistic treatments |],|2|]. Because of indications that small (several percent) variations 
in, for example, neutrino luminosities and shock stagnation radius can make the difference 
in a successful explosion 0, it is clear that even modest relativistic effects comprise an 
indispensible component of realism in supernova studies. 

The multidimensional nature of supernovae must also be recognized as an important as- 
pect of realism. Various observations, especially data from SN 1987A, point to the aspheric- 
ity of supernova explosions (see e.g. || for an overview, and [|J for a recent polarimetry 
analysis of several supernovae). Convection may play an important role in the explosion 
mechanism p|-B|, and the differing results in various simulations in up to two dimensions 
(2D) [P-[T2]J^ JT3|Jl^| show that further study is needed. Ultimately consideration of the third 
dimension will be necessary. Based on an initial exploration of 3D effects, it has been re- 
ported that the sizes of convective cells in 3D simulations are about half as large as in 
2D simulations. Moreover, rapid rotation can significantly affect the strength and spatial 
distribution of convection [Jl5] . Detailed studies of magnetic field generation, jet formation, 
and neutron star kicks also invite 3D treatments. 

Since neutrinos — which carry away about 99% of the gravitational binding energy re- 
leased in the collapse — are believed to drive the explosions, accurate neutrino transport 
is also essential to realistic simulations. Simulations with Boltzmann transport have been 
recently performed in spherical symmetry [ |l6| , |l"7l ,|I|j . 

Including all the physics necessary for realism in a supernova simulation is a daunting 
task. The multidimensional simulations mentioned above, which had simplified neutrino 
transport, taxed the computational resources of their time; the same is true of the recent 
simulations involving Boltzmann transport in spherical symmetry. Adding general relativity 
to the list of desired physics makes things all the more challenging. While numerical relativity 
has been successful in spherical and axisymmetric cases, "...in the general three-dimensional 
(3D) case which is needed for the simulation of realistic astrophysical systems, it has not been 
possible to obtain stable and accurate evolutions...," and it is argued that the difficulties are 



more fundamental than insufficient resolution [T^ . 

In order to overcome the difficulties associated with 3D general relativistic simulations 
and to save resources for 3D hydrodynamics and accurate neutrino transport, an approx- 
imate multidimensional treatment of gravity that captures the phenomena of dynamical 
importance in supernovae would be desirable. The list of new gravitational phenomena in- 
troduced by relativity includes the nonlinearity of the gravitational field, the inclusion of 
all forms of energy and stress as sources, and gravitational waves. The first two of these 
effects are of dynamical importance in supernovae, while gravitational waves will probably 
not exert a strong back-reaction (unless bar or breakup instabilities of some sort become 
operative in the core). The post-Newtonian expansion is systematic and useful in perturba- 



2 



tive applications, but in the present context, an approach that probes the nonlinear nature 
of gravity more deeply would be desirable. A method that could be incorporated into ex- 
isting Newtonian hydrodynamics codes would be even more useful. We here describe such 
an approximation — which we call "NewtonPlus" — and present results for rapidly rotating 
(and therefore highly deformed) neutron stars, which show that this simple "NewtonPlus" 
approach to gravity is indeed a significant improvement over the Newtonian limit. A more 
detailed exposition is given in Ref. JB5 . 



II. EINSTEIN EQUATIONS IN THE NEWTONPLUS APPROXIMATION 

Use of the metric 

ds 2 = -(l + 2$)dt 2 + (1 -2$)dx 2 (1) 

in the Einstein equations yields the Newtonian limit, provided the gravitational potential 
$ < 1 and velocities (including microscopic velocities, large values of which lead to signifi- 
cant stresses) are much less than the speed of light. In order to capture the nonlinearity of 
gravity, the significance of stresses as gravitational sources, and relativistic fluid velocities, 
we propose the use of the following metric: 

where dfl 2 = dO 2 + r 2 sin 2 9d(p 2 . In comparison with Eq. (|I]), the "linearized" metric 
functions have been promoted here to full exponentials, and a second metric function, S, has 
been added. Eq. (^j will then reduce to the Newtonian case if 5 — > 0; we shall see that this is 
in fact the case under conditions in which the Newtonian limit is valid. This provides a tight 
conceptual link with the Newtonian limit. Since it has two independent metric functions, 
this "NewtonPlus" metric should also provide an exact solution in spherical symmetry. 

The metric, a symmetric 4x4 matrix, has ten independent components, but the invari- 
ance of relativity under coordinate transformations implies that in reality there are only six 
degrees of freedom. It is apparent that the NewtonPlus metric of Eq. @ contains only two 
of the six degrees of freedom that should be present. This means that examination of the 
complete set of Einstein equations should reveal inconsistencies, but we have argued that 
these are not too serious in the supernova environment Here we present highlights of 



that more detailed discussion, which involves the (3+1) formulation of the Einstein equations 
PH| in which spacetime is foliated into spacelike slices labeled by a time coordinate. 
We begin with the Hamiltonian constraint. This yields 

V 2 $ = 4ne- 2 *E + l -(d<$>) 2 - L~ m ~ 2i \d t <$>) 2 . (3) 

In this expression V 2 is the usual 3D flat-space Laplacian. The energy density as viewed 
by an "Eulerian" observer (i.e., one whose 4- velocity is orthogonal to the spacelike slices, 
having covariant components = (—a, 0, 0, 0), where a is the lapse function) is denoted by 
E = T^ v n^n v ] where T^ u is the stress-energy tensor. For a perfect fluid, E = T 2 (p + p) — p, 
where T = (1 — v 2 ) -1 ' 2 , v is the magnitude of the fluid velocity as measured by an Eulerian 
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observer, and p and p are respectively the total energy density and pressure in the fluid rest 
frame. We have employed the notation 



OX dY = d r X d r Y + 1 d e X d e Y + 1 «9 X «9 K (4) 

r 2 r 2 sin 

As expected from the conceptual link between the Newtonian and NewtonPlus metrics, Eq. 
(§) identifies $ as a glorified gravitational potential. In addition to the rest energy, the source 
includes internal and kinetic energies and pressure, all boosted by nonlinear contributions 
from $ itself. 

The momentum constraints relate c^<& to the fluid momentum s = T 2 (p + p)v, where 
v is the physical velocity measured by an Eulerian observer. Because the test calculations 
reported here are of stationary configurations, d t § = 0, and the momentum constraints 
don't come into play. The condition d t § = can be taken as a first approximation in the 
supernova environment as well, but it will fail if high density matter moves at relativistic 
velocities. Relativistic velocities may be achieved by infalling matter outside the shock at 
late times, or in outflowing jets or winds; but these situations involve matter at low density 
in comparison with the core. Convection may occur deep in the core, where densities are 
high, but this likely involves nonrelativistic velocities. This means that s = T 2 (p + p)v will 
be arguably small enough everywhere for <9t$ to be neglected. 

Finally we consider the evolution equations of 7^- and Ky. For the NewtonPlus metric, 
explicit calculation shows that the former turn out to be identities leading to no new infor- 
mation. The latter give three different equations for the subdominant metric function 8 that 
are probably inconsistent, though it is not obvious (to us) how to prove this rigorously. It 
turns out that stresses (e.g., pressure) constitute the primary source terms in the equations 
determining S, confirming the expectation expressed previously, that 5 should vanish as the 
Newtonian limit is approached. The observation that 5 will only be appreciable at the high- 
est densities, where pressure begins to make a nontrivial contribution in comparison with 
energy density, suggests a reasonable path forward. In typical cases it is expected that the 
deepest portion of the core will be roughly spherical, even if rapid rotation causes an equa- 
torial bulge of lower density material (e.g., Fig. 2 of Ref. [23]). If this is the case, neglecting 



the angular derivatives of 5 is justified, removing many of the apparent inconsistencies in 
the three equations for 5. The remaining discrepancies have to do with angular derivatives 
of $ and particular components of the stress tensor that appear in each equation. As it 
happens, if one adds all three equations, these remaining discrepancies disappear. Hence 
the equation we shall use to determine 5 is 

d r d r 5 + -d r 5 = 47ce- 2 *{S% + S%) - [«9 r ($ + <5)f , (5) 



r 



where = T h ik hj m , and the spacelike projection tensor is defined by hij = g%j + n^n 
When it is recalled that S e e = S"^ = p in spherical symmetry, this is precisely the equation 
obtained in the spherical case. The source on the right-hand side is to be angle-averaged in 
solving for S. 

It should be straightforward to include the solution of $ and, if desired, 5, in existing 
multidimensional gravitational hydrodynamics codes. The solution of $ would make use of 
the Poisson solver normally used to solve for the Newtonian gravitational potential, the only 
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difference being that one would have to iterate on Eq. @ (with the <9 t $ term dropped) to 
get a self-consistent $. The simplest approximation would be to simply solve for $ in this 
manner, and ignore 5 altogether. The next level of approximation would involve solving Eq. 
(H) for S, but ignoring 5 on the right hand side. Since S is already something of a correction, 
5 appearing on the right hand side is essentially a "correction to the correction." If desired, 
however, Eq. (|5|) could be solved as it stands, with iteration being required. 

The equations of hydrodynamics are obtained from the vanishing divergences of the 
baryon flux vector and stress-energy tensor, employing the NewtonPlus metric. Our study 
of these equations |BJ] shows that they can be cast in a form similar to that used by at least 
one Newtonian multidimensional PPM hydrodynamics code, VH-1. We plan to adapt this 
code to the NewtonPlus approach in the immediate future. 



III. TESTING NEWTONPLUS GRAVITY WITH RAPIDLY ROTATING STARS 

In this section we present calculations of neutron stars undergoing rapid uniform rota- 
tion in order to assess the strengths and weaknesses of the NewtonPlus approximation to 



relativistic gravity. Our models were computed with a code described in Ref. |pq| , which 
was written to compute the structure of relativistic axisymmetric stars. We have modified 
the code to include the ability to perform computations in the Newtonian and various New- 
tonPlus limits: with vanishing metric function 5, with "linearized" 5 (i.e. ignoring 5 on the 
right hand side of Eq. (|5|)), and "full" 5 (solving Eq. (|5|) as written). All of the Newton- 
Plus limits solve a two-dimensional (and stationary) version of the nonlinear Poisson-type 
Eq. (P|) for the enhanced "gravitational potential" $. For the results presented here, the 



high-density portion of the equation of state (EOS) is taken from Ref. pB| |, and is based 



on a field-theoretic description of cold dense matter. We also performed calculations with a 
polytropic EOS of adiabatic index 2, and found qualitatively similar results. 

Panel (a) of Fig. [I] exhibits mass vs. radius curves for spherical stars. It shows that 
while the Newtonian limit exhibits no maximum mass with this EOS,[| the NewtonPlus 
approximation does yield a maximum mass. Even with vanishing 5, the approximation 
captures this consequence of nonlinear gravity. The "linear <5" approximation follows the 
exact relativistic curve until the most dense configurations are reached. Since pressure is the 
main source for 5 (see equation (|J)), the large pressures associated with such high densities 



No turnover in the mass vs. radius curve appears in the Newtonian limit, up to the high-density 
boundary of the tabulated EOS. A configuration with central baryon mass density (total energy 
density) of 3.07 x 10 15 g cm" 3 (4.65 x 10 15 g cm ) has a gravitational mass of 15.6 Mq and 
radius 18.9 km in the Newtonian limit, while the relativistic configuration with this central density 
has a gravitational mass of 1.68 M®. and a radius of 9.40 km. We remind the reader that the 
Chandrasekhar mass phenomenon is a property of stars built on a polytropic equation of state with 
adiabatic index equal to 4/3 (suitable for white dwarfs), but that stars constructed on "realistic" 
nuclear equations of state do not necessarily exhibit this behavior. Instead, the upper mass limit 
of neutron stars derives from the the general relativistic instability indicated by the turning point 
in the mass vs. radius curve. 
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raise 6 to large enough values that it cannot be neglected on the right-hand side of equation 
(H). As expected, the "full 5" approximation is indistinguishable from the relativistic results 
in spherical symmetry, where only two metric functions are needed to describe the spacetime 
exactly. 

Panels (b)-(f) of Fig. [I] show various physical parameters of rapidly rotating configura- 
tions. In order to test the NewtonPlus approximation in a nonspherical setting, we ask the 
question: Given a definite number of baryons rotating at a given uniform angular velocity 
Q, what do the various treatments of gravity do with those baryons? (The fact that baryon 
number is a conserved quantity makes this an obvious way to compare different treatments 
of gravity.) To answer this question we have computed constant baryon mass sequences be- 
ginning at zero rotation (marked by squares) and ending at the mass shedding limit (marked 
by stars). The value of baryon mass chosen, 1.8 M Q , is close to the maximum baryon mass of 
1.95 M Q for the equation of state we employed. The quantities plotted, as a function of the 
(uniform) stellar angular velocity, are gravitational mass; equatorial radius; total angular 
momentum; eccentricity, defined as 1 — r p /r e , where r p and r e are respectively the polar and 
equatorial coordinate radii; and the linear equatorial velocity. 

In panels (b)-(f), the efficacy of the NewtonPlus approximations can be judged by choos- 
ing a value of angular velocity and seeing how close the approximate quantities come to the 
fully relativistic value. While the "full 5" approximation is indistinguishable from full rel- 
ativity in the spherical case, the two curves representing these treatments deviate from 
one another with increasing angular velocity. As expected, the angular velocities at mass 
shedding of the NewtonPlus approximations are closer to the relativistic values than the 
Newtonian case. The NewtonPlus treatments are quite successful at approximating the 
gravitational mass, radius, and eccentricity, while the success of the results for angular mo- 
mentum and equatorial velocity is more modest. (It is expected that the mass and radius 
of the collapsed core are more important to the supernova explosion mechanism than the 
angular momentum.) 

IV. CONCLUSION 

Accurate neutrino transport, 3D hydrodynamics, and relativity are all essential for re- 
alistic supernova simulations. Given the constraints of current hardware, these cannot all 
be treated simultaneously with the detail they deserve. We have therefore presented an 
approximation to full relativity (or set of related approximations) that captures the most 
relevant relativistic effects in the quasispherical supernova environment: nonlinearity creat- 
ing a deeper potential well and pressure being a nontrivial source of gravitation. 

This "NewtonPlus" approach to gravity has a tight conceptual link with the Newtonian 
limit that yields certain advantageous features. The gravitational portion of multidimen- 
sional Newtonian calculations involves only the solution of the Poisson equation for the 
gravitational potential $, and taking its gradient to find the gravitational force. The basic 
idea of our NewtonPlus approach is to promote the Newtonian metric functions — which are 
linear in $ — to full exponentials. We also add a second metric function, 5, whose main 
source is pressure; hence this metric function vanishes in the Newtonian limit. The Ein- 
stein equations yield a nonlinear Poisson-type equation for a (now enhanced) "gravitational 
potential," whose solution in 3D can be obtained in a manner similar to what is currently 



6 



done in the Newtonian limit. The inconsistencies in the Einstein equations arising from the 
reduced number of degrees of freedom turn out to be relegated to the subdominant metric 
function 5; they can be removed by ignoring angular variations in S. This is expected to be 
successful in the supernova context because the region where S makes the greatest difference 
is where pressure is significant in comparison with rest mass density. Normally, this is the 
deepest portion of the collapsed core, which is roughly spherical even when the outer layers 
bulge at the equator due to rapid rotation.^ This strategy — allowing the main contribution 
to the gravitational field to be multidimensional and nonlinear, while allowing a spherical 
correction for the contribution of stresses — reproduces fairly well many of the physical char- 
acteristics of rapidly rotating relativistic stars. Importantly, the hydrodynamics equations 
obtained from the NewtonPlus metric are of the same form as those used in a popular New- 
tonian hydrodynamics algorithm, providing the expectation that existing Newtonian codes 
might be adapted fairly easily to the NewtonPlus approach. 
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2 Ultra-strong magnetic fields |27j or differential rotation (e.g., Ref. |2q| ) can give rise to off-center 
density maxima, making the spherical correction for stresses via the metric function 5 less useful. 
The NewtonPlus approximation with 5 = could still be employed in such (probably exceptional) 
cases, however. 
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FIGURES 



FIG. 1. Panel (a): Mass vs. radius curves for spherical configurations computed with various 
treatments of gravity. Panels (b)-(f): Various physical quantities characterizing uniformly rotating 
configurations, plotted as a function of angular velocity. Each curve represents a constant baryon 
mass sequence computed with the treatments of gravity labeled in panel (a). 
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